
Materials and methods
Data acquisition and preprocessing
All 627 stool metagenomic samples from 11 published studies were downloaded from the NCBI/EBI databases using Kingfisher v0.4.1 [Woodcroft, 2022]. Data quality was assessed with FastQC v0.12.1 [Andrews et al., 2023], and raw reads were processed with fastp v0.23.4 to remove low-quality sequences [Chen et al., 2018]. Human-derived sequences were filtered out using HISAT2 v2.2 [Kim et al., 2019] with the GRCh37 human genome (Release 47) as a reference [Schneider et al, 2017].
Assembly and binning
Metagenomic assembly was performed using MEGAHIT v1.2.9 [Li et al., 2015], retaining contigs longer than 1000 bp. Filtered contigs were aligned to metagenomic reads using HISAT2 v2.2 [Kim et al., 2019]. Binning was conducted in two stages: (1) initial binning with MetaBAT2 v2.12.1 [Kang et al., 2019], MaxBin2 v2.2.7 [Wu et al., 2016], and Semibin2 v2.1.0 [Pan et al., 2023], followed by (2) refinement with DAS Tool v1.1.7 to improve bin quality [Sieber et al., 2018]. Bins were quality-checked with CheckM v1.2.3 [Parks et al., 2015] and dereplicated at 98% nucleotide identity using dRep v3.4.5 [Olm et al., 2017] to generate operational genomic units (OGUs).
Taxonomic profiling and analysis
Taxonomic annotation of bins was performed using GTDB-Tk v2.1.1 [Chaumeil et al., 2022] with the GTDB r207 database [Parks et al., 2022]. Additional data processing utilized Samtools v1.17 [Li et al., 2009], BEDTools v2.31.0 [Quinlan et al., 2010], and BBMap v39.06 [Bushnell et al., 2014]. Finally, OGU abundance profiles were generated by mapping reads to the OGU catalog using HISAT2 v2.2, followed processong with InStrain v1.9.0 [Olm et al., 2021]
Marker OGU discovery
The identification of OGUs associated with immunotherapy outcomes followed an established analytical framework from our previous works [Olekhnovich et al., 2023; Zakharevich et al., 2024]. Differential rankings was first performed using Songbird [Morton et al., 2019] to identify OGUs showing relative abundance variations between experimental groups, applying a conservative absolute differential value threshold of > 0.3. For candidate OGUs meeting this criterion, we subsequently calculated log-ratio abundances using Qurro [Fedarko et al., 2020] and determined statistical significance through Wilcoxon rank-sum tests implemented in the R statistical environment. The biomarker selection process incorporated stringent cross-validation criteria to ensure robust identification of microbial signatures. OGUs demonstrating consistent positive associations with therapeutic response across multiple datasets were retained as potential beneficial biomarkers, while any evidence of negative association with treatment outcome in any dataset resulted in automatic exclusion regardless of other positive associations.This approach enabled simultaneous identification of two clinically meaningful biomarker categories: microbial taxa positively correlated with successful immunotherapy outcomes and those associated with adverse therapeutic responses. The methodology emphasizes reproducibility through multi-dataset validation and maintains rigorous standards for biomarker qualification by requiring consistent directional effects across independent cohorts.
Markers testing
Determined R/NR biomarkers were used to calculate log ratios, followed by statistical assessment using the lmer function from the lmerTest package [Kuznetsova et al., 2017] in the R statistical environment. Post-hoc testing was performed using Tukey’s method with the multcomp package [Bretz et al., 2011], applying Bonferroni correction for multiple comparisons. Model residuals were tested for normality using the Shapiro-Wilk test (shapiro.test function in R). Additional testing of marker OGU included investigating the ability of calculated log ratios to separate experimental groups using logistic regression using the method of leave-one group out cross validation. Gene-set enrichment analysis (GSEA) implemented in GSEA function from clusterProfiler package [Xu et al., 2024] using to investigate statistical relationships of OGU marker sets across taxonomy and other attributes included food, body site or pathogens.GSEA analysis results were visualized using gseavis package [Zhang et al., 2025].
Functional profiling of marker OGU
Functional profiles of OGUs were generated using MetaCerberus [Figueroa et al., 2024] with annotations from the Carbohydrate-Active EnZymes (CAZy) database [Drula et al., 2022] and the Kyoto Encyclopedia of Genes and Genomes (KEGG) [Kanehisa et al., 2025]. Statistical associations between functional features and experimental groups were assessed using logistic regression. Gene set enrichment analysis (GSEA) was applied to evaluate links between gene sets and immunotherapy response, with results visualized using the ggseavis package.
Additional analysis and visualization
For visualization and statistical analysis, we utilized the following R packages: ggplot2 [Wickham et al., 2016] for graphics creation, ggpubr [Kassambara, 2023] for publication-ready plots, hrbrthemes [Rudis, 2020] for enhanced themes, and RColorBrewer [Neuwirth, 2022] for color palette management. String data manipulation was performed using the stringr package [Wickham, 2022]. The data analysis report was generated with Quarto [Allaire et al., 2022], incorporating dynamic tables via DT [Xie et al., 2023] and enhanced download functionality using downloadthis [Sidi, 2023].
Data availability
In this study, we used open access data from the NCBI/EBI Sequence Read Archives, identified by the following BioProjects accession numbers: PRJNA397906, PRJEB22893, PRJNA399742, PRJNA770295, PRJEB43119, PRJNA762360, PRJNA1011235, PRJNA615114, PRJNA866654, PRJNA494824, PRJEB49516. MAGs catalog assembly pipeline were description in https://github.com/JeniaOle13/Cancer_MAGs. All MAGs sequences were deposited in NCBI under accession PRJNA1196825. Songbird Qurro files for QIIME2 viewer were available at Zenodo. Source code and Quarto report available at https://github.com/JeniaOle13/cancer-biomarkers.
References
- Chen, Shifu, et al. “fastp: an ultra-fast all-in-one FASTQ preprocessor.” Bioinformatics 34.17 (2018): i884-i890.
- Li, Dinghua, et al. “MEGAHIT: an ultra-fast single-node solution for large and complex metagenomics assembly via succinct de Bruijn graph.” Bioinformatics 31.10 (2015): 1674-1676.
- Kim, Daehwan, et al. “Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype.” Nature biotechnology 37.8 (2019): 907-915.
- Kang, Dongwan D., et al. “MetaBAT 2: an adaptive binning algorithm for robust and efficient genome reconstruction from metagenome assemblies.” PeerJ 7 (2019): e7359.
- Wu, Yu-Wei, Blake A. Simmons, and Steven W. Singer. “MaxBin 2.0: an automated binning algorithm to recover genomes from multiple metagenomic datasets.” Bioinformatics 32.4 (2016): 605-607.
- Pan, Shaojun, Xing-Ming Zhao, and Luis Pedro Coelho. “SemiBin2: self-supervised contrastive learning leads to better MAGs for short-and long-read sequencing.” Bioinformatics 39.Supplement_1 (2023): i21-i29.
- Sieber, Christian MK, et al. “Recovery of genomes from metagenomes via a dereplication, aggregation and scoring strategy.” Nature microbiology 3.7 (2018): 836-843.
- Parks, Donovan H., et al. “CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes.” Genome research 25.7 (2015): 1043-1055.
- Olm, Matthew R., et al. “dRep: a tool for fast and accurate genomic comparisons that enables improved genome recovery from metagenomes through de-replication.” The ISME journal 11.12 (2017): 2864-2868.
- Chaumeil, Pierre-Alain, et al. “GTDB-Tk v2: memory friendly classification with the genome taxonomy database.” Bioinformatics 38.23 (2022): 5315-5316.
- Parks, Donovan H., et al. “GTDB: an ongoing census of bacterial and archaeal diversity through a phylogenetically consistent, rank normalized and complete genome-based taxonomy.” Nucleic acids research 50.D1 (2022): D785-D794.
- Li, Heng, et al. “The sequence alignment/map format and SAMtools.” bioinformatics 25.16 (2009): 2078-2079.
- Quinlan, Aaron R., and Ira M. Hall. “BEDTools: a flexible suite of utilities for comparing genomic features.” Bioinformatics 26.6 (2010): 841-842.
- Bushnell, Brian. “BBMap: a fast, accurate, splice-aware aligner.” (2014).
- Olm, Matthew R., et al. “inStrain profiles population microdiversity from metagenomic data and sensitively detects shared microbial strains.” Nature Biotechnology 39.6 (2021): 727-736.
- Olekhnovich, Evgenii I., et al. “Consistent stool metagenomic biomarkers associated with the response to melanoma immunotherapy.” mSystems 8.2 (2023): e01023-22.
- Zakharevich, Natalia V., et al. “Systemic metabolic depletion of gut microbiome undermines responsiveness to melanoma immunotherapy.” Life Science Alliance 7.5 (2024).
- Morton, James T., et al. “Establishing microbial composition measurement standards with reference frames.” Nature communications 10.1 (2019): 2719.
- Fedarko, Marcus W., et al. “Visualizing’omic feature rankings and log-ratios using Qurro.” NAR genomics and bioinformatics 2.2 (2020): lqaa023.
- Bolyen, Evan, et al. “Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2.” Nature biotechnology 37.8 (2019): 852-857.
- Kuznetsova, Alexandra, Per B. Brockhoff, and Rune HB Christensen. “lmerTest package: tests in linear mixed effects models.” Journal of statistical software 82 (2017): 1-26.
- Bretz, F., Hothorn, T., & Westfall, P. (2011). Multiple comparisons using R. CRC Press.
- Xu S, Hu E, Cai Y, Xie Z, Luo X, Zhan L, Tang W, Wang Q, Liu B, Wang R, Xie W, Wu T, Xie L, Yu G (2024). “Using clusterProfiler to characterize multiomics data.” Nature Protocols, 19(11), 3292-3320. doi:10.1038/s41596-024-01020-z
- Zhang, Jun, et al. “GseaVis: An R Package for Enhanced Visualization of Gene Set Enrichment Analysis in Biomedicine.” Med Research (2025).
- Figueroa III, Jose L., et al. “MetaCerberus: distributed highly parallelized HMM-based processing for robust functional annotation across the tree of life.” Bioinformatics 40.3 (2024): btae119.
- Drula, Elodie, et al. “The carbohydrate-active enzyme database: functions and literature.” Nucleic acids research 50.D1 (2022): D571-D577.
- Kanehisa, Minoru, et al. “KEGG: biological systems database as a model of the real world.” Nucleic Acids Research 53.D1 (2025): D672-D677.
- Wickham, H. (2016). ggplot2: Elegant graphics for data analysis (2nd ed.). Springer-Verlag.
- Kassambara A (2023). ggpubr: ‘ggplot2’ Based Publication Ready Plots. R package version 0.6.0.
- Rudis, Bob. hrbrthemes: Additional Themes, Theme Components and Utilities for ‘ggplot2’. R package version 0.8.0, 2020.
- Neuwirth, Erich. RColorBrewer: ColorBrewer Palettes. R package version 1.1-3, 2022.
- Wickham, Hadley. stringr: Simple, Consistent Wrappers for Common String Operations. R package version 1.5.0, 2022.
- Allaire, J. J., et al. Quarto. Version 1.3, 2022.
- Xie, Yihui, Joe Cheng, and Xianying Tan. DT: A Wrapper of the JavaScript Library ‘DataTables’. R package version 0.28, 2023
- Sidi, Fahim. downloadthis: Implement Download Buttons in ‘R Markdown’. R package version 1.0.0, 2023
- Woodcroft, Ben J. Kingfisher: Rapid SRA Download Tool. Version 0.4.1, 2022
- Andrews S. FastQC: A quality control tool for high throughput sequence data. Version 0.12.1, 2023.
- Schneider VA, et al. Evaluation of GRCh38 and de novo haploid genome assemblies demonstrates the enduring quality of the reference assembly. Genome Research 2017;27(5):849-864.
Results
Data overview
The analysis incorporated 627 gut metagenomic profiling samples obtained from 11 independent external datasets. Patients were stratified by immunotherapy response into two groups: responders (R group, n = 365; 58.2%) and non-responders (NR group, n = 262; 41.8%). Response assessment followed RECIST 1.1 criteria, with the R group including patients showing complete response (CR), partial response (PR), or stable disease (SD) at 6-month follow-up, while the NR group comprised exclusively progressive disease (PD) cases. The study cohort received various immunotherapy regimens, including anti-PD1, anti-CTLA4, or combination therapies. Cancer type distribution revealed melanoma predominance (n = 456; 72.7%), followed by gastrointestinal (GI) cancers (n = 82; 13.1%), non-small cell lung cancer (n = 15; 2.4%), breast cancer (n = 4; 0.6%), ovarian cancer (n = 2; 0.3%), and other malignancies (n = 68; 10.8%). All samples were collected prior to treatment initiation to evaluate baseline microbiota status.
Quality control
Following quality filtering with fastp, aggregated quality metrics were compiled using MultiQC to summarize the FastQC reports (MultiQC are attached to the report below). Before filtering, approximately ~15bn reads (mean ± SD: 24 ± 15 mln reads per sample) were retained for downstream analysis, with an average of 34 ± 18% of reads per sample removed during quality control (see Table 2 for detailed filtering statistics).
Taxonomic composition
A total of 3649 OGUs belonging to 12 bacterial phyla were detected in 627 samples. The most numerous of these were Firmicutes (n = 2443), Actinobacteriota (n = 590), Bacteroidota (n = 394), Proteobacteria (n = 140), Desulfobacterota (n = 27) and others (n = 55). OGU abundance profiles of stool samples presented in Table 3.
Marker OGUs discovery
The analysis of marker OGUs using Songbird and Qurro identified 793 significant features, with 326 enriched in responders (R) and 467 enriched in non-responders (NR). Among R markers, the most prominent species included Blautia wexlerae (n = 49), Faecalibacterium prausnitzii (n = 27), Fusicatenibacter saccharivorans (n = 15), Gemmiger qucibialis (n = 14), and Blautia faecis (n = 8). In contrast, non-responders exhibited enrichment of different microbial taxa, including F. prausnitzii (n = 17), Dysosmobacter sp001916835 (n = 12), Dialister invisus (n = 11), Veillonella parvula (n = 9), and ER4 sp000765235 (n = 7). Notably, while F. prausnitzii was detected in both groups, its higher abundance in responders (27 vs. 17 OGEs) implies that strain-level differences or functional variations could influence clinical outcomes. The strong association of Blautia species, particularly B. wexlerae, with responders highlights its potential as a biomarker for positive treatment response. It is also important to note the association of oral bacteria V. parvula and D. invisus with negative outcomes of immunotherapy. Marker OGEs description as detailed in Table 4.
Saving 8 x 3 in image


lmer test results
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: lmer(formula = log_ratio ~ response + (1 | dataset), data = df.log_ratio)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
R - NR == 0 1.3472 0.1286 10.48 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- bonferroni method)
Dataset variable inportance
ANOVA-like table for random-effects: Single term deletions
Model:
log_ratio ~ response + (1 | dataset)
npar logLik AIC LRT Df Pr(>Chisq)
<none> 4 -979.08 1966.2
(1 | dataset) 3 -986.73 1979.5 15.286 1 9.241e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residuals of the model normality testing
Shapiro-Wilk normality test
data: residuals(lmer.log_ratio)
W = 0.99234, p-value = 0.007271
Effect size estimation
Cohen's d | 95% CI
--------------------------
-0.89 | [-1.07, -0.71]
- Estimated using pooled SD.


Marker OGUs features



